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Abstract 

The outcome of competition among species is influenced by the spatial distribution 
of species and effects such as demographic stochasticity, immigration fluxes, and the 
existence of preferred habitats. We introduce an individual-based model describing 
the competition of two species and incorporating all the above ingredients. We find 
that the presence of habitat preference — generating spatial niches — strongly sta- 
bilizes the coexistence of the two species. Eliminating habitat preference — neutral 
dynamics — the model generates patterns, such as distribution of population sizes, 
practically identical to those obtained in the presence of habitat preference, pro- 
vided an higher immigration rate is considered. Notwithstanding the similarity in 
the population distribution, we show that invasibility properties depend on habitat 
preference in a non-trivial way. In particular, the neutral model results results more 
invasible or less invasible depending on whether the comparison is made at equal 
immigration rate or at equal distribution of population size, respectively. We discuss 
the relevance of these results for the interpretation of invasibility experiments and 
the species occupancy of preferred habitats. 

Key words: Spatial models. Dispersal, Voter model. Heterogeneous habitat, 
Neutral Theory 



1 Introduction 



A central problem in community ecology is to understand the ecological forces 
leading to the observed patterns of coexistence or exclusion of competing 
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species ( jRicklefs and Schluted . Il993l : iBrown et al.l . Il995l ). This issue is im- 
portant f o r und erstanding: both simple communities made up of few species 
( IChessoru . 120001) and "biodive rsity hotspots" with a large number of coexist- 
ing species ( iLeigh et al.l . 120041 ) . Historically, this problem has been approached 
at two distinct levels. On the one hand, focus has been put on the detailed 
mechanisms of interaction between species (e.g. intra- and inter-specific com- 
petitions) caused by their differenti ation in exploiting resourc es, resulting in 
the concept of the ecological niche ( IChase and Leiboldl. 120031). For example, 
it has been shown how habitat heterog eneity (jBeckage and Clarkl . |2003[ ) or 
a tradeoff in dispersal range strategies ( iBolker and Pacalal . Il999[ ) may pro- 
mote coexistence. An alternative explanation for the observed species richness 
and distribution is in terms of processes intrinsically due to chance, su ch as 
colonization, immigration and extinction (IMacArthur and Wilson! . 119671 ). dis- 
regarding differences among species. 



In recen t years , the neutral theory of biodiversity (iHubbelll . Il979l : iBelll . 12001 
Hubbelll . l200l[ ) considerably developed the latter approach by explicitly as- 
suming equivalence among species at the individual level. The interest in the 
neutral theory has been triggered by its ability to successfully predict several 
biodiversity patterns observed in tropical forests, s uch as species-abundance 



distributions in diff erent permanent sampling p lots (IBelll. l200ll: iHubbelll. 12001 



Volkov et al.l. 120031) a nd species-area relations (IDurrett and Levinl . Il996t iBell 



200ll : IHubbelll . 1200 ll ). Its success underlined the importance of stochasticity 
(ecologica l drift) and dispersal limitation in the assemblage of natural com- 
munities ( Chave . 2004 : Alonso et al. . 20061). wh i ch are now recognized as key 



ingredient also in niche-based models ( iTilmanl . |2004| ). However, niche-based 



to neutral ones when compared with data (I Chave et al 



Mouquet and Loreaul . l2003l : iGilbert and Lechowiczl . 12004 
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). This 



suggests that these patterns tend to average out the dependence on the de- 
tails of the theory (see also the discussion in Pueyo et al. ( 2007f )) and thus 



cannot be used to discriminate the relative importance of niche-based and 
neutral forces. In this perspective, the study of dynamical prope rties such as 



invas ibility can be a promising way to disentangle these effects (iDaleo et al. 
20091 ). 



To understand the key differences between neutral and non-neutral compe- 
tition, it is useful to consider models that can be continuousl y tuned from 



niche - based to neutral settings by v a rying some parameters (IChave et al 



20021 : iGravel et al.l . 120061 : lAdler et al.l . l2007l ). An obvious difficulty with this 



pro gram comes fr om the unavoidable complexity of realistic niche-based mod- 
els ( IChasd . I2OO5I ). where species are not equivalent and the environment is 
heterogeneous both in space and time. This suggests an approach whereby 
simplified models with few parameters are studied, for example by making 
some specific assumptions on how neutrality is violated. 



2 



In this paper we study the dynamics of two species A and B that compete for 
space. The model is devised in such a way that a single parameter controls the 
overlap between the niches occupied by the two species, from complete - neu- 
tral - to no overlap - two independent niches. The model is individual-based 
and incorporates the basic ingredients of neutral theory: coexistence results 
from immigration from a metacommunity, balancing demographic stochastic- 
ity, which alone would lead to extinction. Niches are introduced in this neutral 
scenario via preferential habitats: half of the sites in the ecosystem are favor- 
able for the colonization of individuals of one species and the other half are 
favorable for the other species. The ecologica l advantage is realized through a 



biased "lottery" (IChesson and Warnerl . Il98ll ). We consider a symmetric situ- 



ation by choosing the same statistical bias, 7, for individuals of species A and 
B to colonize their respective preferred habitats. When 7 = (no habi t at di- 



versification) the rnodel reduces to the voter model (IHoUey and Liggettl . Il975 



Cox and Griffeathl . ll986l ). which is a prototype of neutral dynamics. Increasing 
7, species acquire an advantage in colonizing some sites and a disadvantage 
in others. A very large 7 eventually leads to segregation of the two species to 
their preferential habitats. Segregation will be complete when the choice of 
dispersal allows individuals to reach all their preference sites or incomplete in 
the presence of dispersal limitation. 

The aim of this work is to use this simple model to understand the effect 
of habitat diversification on coexistence and dynamics of ecological commu- 
nities. In particular, our concern will be on contrasting the effect of habitat 
diversification with the neutral model where no preferred habitat exists. 



2 Model 



We consider an individual based, spatially-explicit model of a community made 
of two competing species A and B with population A^^ and Nb, respectively. 
The community lives in a patch made oi N = sites on a square lattice of 
side L, on which we assume periodic boundary conditions. Each lattice site 
is occupied by a single individual of one of the two species. For the sake of 
simplicity, we assume that the patch is saturated, i.e. with no empty sites 
— each dead individual is immediately replaced, so that the total number 
of individual is constant, N a + -/Vr = A^. The latter hypothesis is commonly 



assumed for its convenience (IHubbelll . I2OOII : IChave et al.l . |2002| ) and, strictly 
speaking, corresponds to considering infinite fecundity. However, a finite but 
reasonably high fecundity wou ld lead to almost-satura t ed ecosystems wit h 
qualitatively similar dynamics ( Durrett and Levin . 1996 : Chave et al.l . I2OO2I ). 



Our main interest here is to study the effect of habitat preference on compe- 
tition. To this aim, we assign to each site a specificity: half of the sites are 
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favorable (as below specified) to individuals of species A and the other half to 
individuals of species B. We denote such sites by a and b, respectively. The 
site specificity can have several different (often concomitant) ecological ori- 
gin s such as abundance of resources, predation pressure (see, e.g., the review 
by lAmarasekard . |2003[ ). and/or environmental con ditions such as elevatioii , 
temperature, soil moisture or other parameters as inlZillio and ConditI (120071 ) 
and as suggested by observations (IBeckage and Clarkl . |2003| ). The net effect 
of these different mechanisms is here assumed to increase by a factor 7 the 
chance of individuals to colonize a preferred site. This is illustrated in the top 
panel of FigJH to be compared with the bottom cartoon which shows the neu- 
tral model, without site specificity. Simi lar models have been pro posed also 
in the cont ext of heterogeneous catalysis (IFrachebourg et al.l . 119951 ) and social 
dynamics (IMasuda et al.l . I2OIOI ) . 



For the sake of simplicity, site specificity is randomly assigned at the begin- 
ning and left unchanged during the dynamics. Of course, in natural ecosys- 
tems, spatial arrangement of sites with a certain specificity will usually be 
characterized by a certain degree of correlations, which will in general tend 
to enhance niche effects. In this respect, we expect that our choice will tend 
to underestimate the effect of habitat preference. Clearly, the model can be 
generalized by introducing asymmetries, i.e. different 7's for the two species 
or different fractions of advantageous sites. Here we shall limit our analysis 
to the simple symmetric case, so that no species has a net advantage and the 
degree of habitat preference is controlled by a unique parameter. 

We also assume a continuous immigration in the patch of individuals A ot B 
at rate i'. This inflow is necessary to avoid the drift to extinction of one of the 
two species. 

For any given size L of the patch, which flxes the number of individuals = 
L^, the model is controlled by two parameters only: the colonization advantage 
7 and the immigration rate u. The elementary step of the dynamics is as follows 

i) a site is randomly chosen and the individual there residing is killed; 

ii) with probability (1 — z/), the individual is replaced by a copy of one of the 
four neighbors, chosen via a lottery which gives a competitive advantage 
(modeled as a weight 7) to individuals having that site as preferred habitat 
(see the sketch in Fig{T]and eqs. ([1]) and ([2]) below). 

iii) with probability z/, the individual is replaced by an immigrant. For sim- 
plicity, we assume the two species being equipopulated at the metacom- 
munity level, so that the probability of being replaced by an individual of 
one of the two species is 1/2, apart from the competitive advantage on the 
speciflc empty site. 

In formulas, steps ii) and iii) can be expressed as follows. If the individuals 
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Fig. 1. Sketch illustrating the model with (top) and without (bottom) habitat pref- 
erence. Two representative configurations of 4 x 4 systems are shown, white squares 
are advantageous to A, gray ones to B, and on the right we sketch the lottery 
dynamics (eqs. ([1] and ([2])): the width of arrows represents the habitat preference 
intensity 7. Notice that in the bottom panel the width of arrows is insensitive to 
the site specificity 

is killed in a site of type a advantageous for individuals of species A, the 
probabilities of being replaced by an individual A or B are given by 



'I 



(1+7)"A _|_ j, , 1+7 
(l+7)nA+ns 2+7 



W%{nA.nB) 



m L 

(l+7)nA+nfl ^ 2+7 



respectively, where ua and ub denote the number of individuals of species A 
and B in the neighborhood of the considered site. Similarly, if the individual 
is killed in a site of type 6, we have 



W\{nA,nB) 



'1 



nA + (l+7)"s 



+ Z/ 



2+7 



Wl,{nA,nB) 



I (1+7)"B _|_ , , 1+7 
n..4 + (l+7)"s 2+7 



The competitive lottery ([T]) and ([2]) used in step ii) r epresents a bia sed (non- 
neutral) generalization of that used in neutral models ( lHubbelll . l200ll ). Indeed, 
for 7 = 0, the neutral dynani i cs of the voter model with only two species is 



recovered (IHolley and Liggettl . Il975l ) . 



Notice that the model is set up in such a way that the fitness advantage to be 
on colonization: after having colonized a site, mortality and dispersal do not 
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depend on being on a preference site. In other terms, the fitness advantage 
belongs to the seeds and not to the individuals th e msely es. In this interpre- 
tation, and at variance with IChesson and Warner! (119811 ). we excluded from 
the contribution to the (implicit) seed pool of the individual who died. Al- 
though th e latte r may be more realistic in modeling, e.g., perennial plants 



( iLin et al.l . |2009| ) , we preferred the first to directly compare the ne utral ver- 
sion of the model (7 = 0) originally proposed by iHubbelll (119791 ) with the 
non-neutral (7 > 0) variants. We do not expect big differences in the out- 
come of the two models, a part from a weak discrepancy in the percentage of 
occupancy of preferential habitats. We also note that this difference becomes 
less relevant when longer dispersal is introduced. Indeed, the model can be 
generalized to different form of dispersal mechanism s, e.g. with a given fi- 
nite r ange. Previous investiga t ions on similar models (IRosindell and Cornell! . 
2OO7I : IZillio and Conditl . boOTi : IPigolotti and Cencinil . l2009|) have shown that 
different dispersal mechanisms yield qualitatively similar results, as far as they 
are finite-ranged and all species adopt the same dispersal strategy. Here, we 
consider the two limiting situations of nearest-neighbor and global dispersal. 



In the global dispersal case, all individuals present in the patch can colonize 
the site left empty by the dead one. This limiting case strongly simplifies the 
simulations of the model since it makes the distance and spatial distribution 
of the sites irre levant. The gl obal dispersal model may be thought as a two 
islands model (jWrightl . Il93ll ). where each of the two islands has room for 
N/2 individuals and contains all the sites favorable to one of the two species. 
The state of the system is then unequivocally identified by the numbers NAa 
and NBh of individuals of the two species occupying the respective sites of 
preference. All the other variables can be expressed in terms of N^a and Nsb-, 
for instance: the number of individuals of species A (resp. B) outside their 
preference sites is NAb = N/2 — Nbb (resp. Nsa = N/2 — NAa) and the total 
number of individuals of species A (resp. B) is Na = N/2 + NAa — Nbb (resp. 
A^^ = N/2 + NBh — ^Aa)- The evolution of the system is thus determined by 
the probabilities per elementary step that NAa and NBb increase or decrease 
by one unity: 



WNBb^NBb + l 



\-^) WI{Na.Nb) 



Bb 



1 _ A^ 

2 A^ 

NAa 

N 

NBb 

N 



W^Na^Nb) 



W%{Na,Nb) 
W\{Na,Nb). 



(2) 



where Wy is given by eqs. ([T]) and (j2]) with ua and hb replaced by A^i and 
Nb-i respectively. 
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For 7 = 0, the model is neutral and reduces to the ordinary voter model with 
on ly two species and i mmig ration (the multi-species version was considered 



m 



Durrett and LevinI (119961)) in th e spatially explicit case, or to the Moran 



model with mutation (Moran, 19581 ) in the global dispersal version. Conversely, 
when 7 becomes very large, the competitive advantage gets so intense that 
individuals tend to localize in their preferred habitats, with an extremely low 
chance of colonizing the rest of the ecosystem, so that the two species do not 
compete anymore. 



3 Results 



Extinction times and fixation probabilities 



In the absence of immigration (z/ = 0) and for any choice of A^, 7 and dispersal 
range, persistent coexistence is not possible, since species cannot recover from 
a local extinction event. Ecological drift will ultimately drive one of the two 
species to extinction, and the dominance of the other will become stable — 
reachi ng fixation as from population genetics terminology, see e.g. I Gillespie 
( 119941 ). Studying the propert ies of the dynamics toward e xtinction i s how ever 



interesting and informative (IChesson and Warnerl . Il981t IChesson I . Il982l ). In 



particular, the time needed for the extinction of one of the species is important 
to understand how crucially biodiversity depends on a steady immigration flux. 
Moreover, the way the extinction time and the probability of fixation of one 
of the two species depend on the deviation from neutrality allow to quantify 
to what extent habitat preference promotes coexistence. 



We start describing a set of simulations in which the two species are equipop- 
ulated at the initial time, i.e. A^^ = A'"^ = N/2, and individuals are randomly 
placed in the system. The dynamics is then followed until the extinction of 
one of the species. By averaging over many realizations of the dynamics (from 
5 X 10^ to 10^), we can access the extinction time. Here we focus on the 
average time, T, tho ugh also the fluctuations have an ecological relevance 
(jPigolotti et all . |2005[ ). The average extinction time T as a function of the 
community size N and for several values of 7 is shown in Figj2] as obtained 
with the global dispersal model. As customary in models with overlapping 
generations, we measure T in generations, where the time unit corresponds to 
A^ successive time steps i)-iii) described in the m odel section. F or 7 = 0, we 
recover the result expected for the Moran model (jMoraru . 119581 ) i.e. T (x N. 
While the presence of habitat preference 7 > leads to a dramatic increase 
of the average extinction time: for large enough populations, T becomes ex- 
ponentially large in the community size A^, i.e. 
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Fig. 2. Results for the model with global dispersal without immigration (v = 0): (a) 
Mean extinction times T as a function of N for different 7, as in label. The black 
lines indicate exponential fits of the form T = C exp{C'{'y)N) while the gray straight 
line indicate the neutral expectation T (x N. The inset shows C'{'y) versus 7 for 
both the global dispersal (GD) and nearest neighbor (NN) models, (b) Probability 
of fixation Pfix of A vs the initial fraction x = Na/N for = 40. Inset: Pfix vs 
holding fixed x = Na/N = 0.1. 

where C depends on 7 (as shown in the inset of Figl2|L). The exponential 
dep endence implies that the coexistence may be considered stable for large 



(jChesson I . Il982l ). The spatially explicit version of the model (not shown 



here) presents the same qualitative features with small quantitative changes: 
i n the neut r al cas e 7 = 0, logarithmic corrections are present T oc A^log(A^) 
(iKrapivskyl . Il992l ) and the exponential rates for 7 > are slightly different 
from those obtained with global dispersal (inset of Figj2|i). In large communi- 
ties, even a tiny habitat preference leads to a dramatic increase in the average 
extinction time, stabilizing the system on any realistic timescale. 



To further confirm the stabilizing effect of habitat preference, we consider a 
different setting in which species A is present at initial time with a fraction 
of individuals x = N^/N. By averaging over 10^ realizations for each x, we 
computed the probability that sp ecies A become s fixated, Pfix{x) (Fig|2b). 
In the neutral case, Pfix{x) = x ( jGillespid . 119941 ) . while Pfix{x) is closer to 
1/2 when 7 is increased, the effect being stronger the larger is the community 
size (see inset). Therefore, strong habitat preference tends to compensate any 
initial disproportion between two large populations, making equally likely the 
fixation of one of the two species. In particular, if A is the invading species, i.e. 
the less represented at the beginning (x <^ 1/2), its chances of invading the 



systems greatly increase in the pres ence of habitat preferen ce, in agreement 
with other models and observations (IMelbourne et al.l . 120071 ). 



Regimes of coexistence 



The presence of an immigration pressure allows species to recover from local 
extinction events, ensuring dynamical coexistence, as illustrated in the left 
panels of Fig. [3l where the evolution of the population Na is shown over 
10^ generations at varying the rate of immigration u. The figures refer to 
the global dispersal model with a habitat-preference intensity 7 = 0.3. We 
mention that qualitatively similar features are observed also in the nearest- 
neighbor dispersal case and also when the model is neutral. 



As expected from classical theories ( iMacArthur and Wilson! . Il967( ). increas- 
ing the immigration pressure enhances the degree of coexistence of the two 
species. At fixed N, if u is small, one of the two species dominates for most 
of the time. Stochastic fluctuations, whose amplitude decreases at increasing 
N, lead to an alternation in the dominating species (top panel of Fig. E]). The 
turnover in the dominating species takes place on timescales rapidly growing 
with N (roughly of the order of the average extinction time). At intermedi- 
ate values of z/, coexistence is possible, but episodic local extinctions can rule 
out a species from the system for several generations (middle panel). Finally, 
increasing u even further, local extinctions cannot persist and coexistence be- 
comes the rule: it is more and more likely to observe states in which the two 
species are roughly equipopulated (bottom panel). We quantify these three 
behaviors in terms of the functional shape of the probability P{Na), averaged 
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Fig. 3. Different regimes of coexistence in the presence of immigration in the global 
dispersal model with habitat preference: evolution of the population Na across 10'^ 
generation for a system of = 100 individuals (left) and the corresponding distri- 
butions P{Na) (right). The three panels are obtained holding the habitat-prefer- 
ence intensity fixed 7 = 0.3 and varying the immigration rate v. (top) u = 0.003 
corresponding to monodominance; v = 0.015 leading to a mixed phase; ly = 0.05 
displaying coexistence (see text for details). 
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Fig. 4. Phases in the N — u parameter space. White regions denote monodominance, 
gray-filled regions pure coexistence, and dashed regions mixed phases. On the left 
(panels a-c) data obtained with the global dispersal model, while on the right (panels 
d-f ) with the nearest-neighbor model. The three rows correspond to different values 
of 7 = 0, 0.3, 1.0 (from top to bottom). Axis are in log scale. Notice that for 7 = 1 
and global dispersal (panel c) we could not show the curves for > 240 due to lack 
of statistics. 



over time, of one of the two species to oc cupy a given fraction of the ecosys- 
tem (see, e.g. iLoreau and Mouquetl . Il999l ). We recall that there is a complete 
symmetry between A and B so that P{Na) = P{N — Na) = P{Nb)- In par- 
ticular, with reference to the right panels of Fig. [3], from top to bottom we 
can distinguish the following classes of distributions corresponding to different 
regimes of coexistence, as determined by the number of maxima of P{Na)'- 
Monodominance: P{Na) is U-shaped, achieving its two maxima at Na = 
and Na = N; 

Mixed: P{Na) has three maxima at Na = and Na = N a.s before plus 
Na = N/2, meaning alternation between states of monodominance and pure 
coexistence as defined below; 

Pure coexistence: P{Na) is bell-shaped and has a single maximum at Na = 
N/2, so that most of the time the populations fluctuate around the equipop- 
ulated state. 



We quantify the stabilizing effect of habitat preference by studying for which 
values of the parameters 7, z/, and N, the above defined regimes are observed. 
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Figure H] shows, for different values of 7, how the monodominance (white), 
coexistence (gray) and mixed (patterned) regimes organize in the N — u plane 
for both global (left) and nearest-neighbor (right) dispersal. 

As observed for the extinction times, nearest-neighbor and global dispersal 
models display qualitatively similar features. The only difference is that in the 
neutral case 7 = the global dispersal model never shows a mixed distribution; 
conversely, with nearest-neighbor dispersal, a tiny region of mixed phase exists 
also for 7 = 0, compare Fig. and d. For 7 = and global dispersal, 
the transition between monodominance and coexistence can be analytically 
obtained and occurs for z/ > Uc{N) = 2/(2 + A^) (see Appendix |A|) . 

The most surprising feature is that the critical immigration rate vdN) for 
observing coexistence seems to be independent of 7 and of the dispersal. In- 
deed, we could not determine appreciable quantitative differences between the 
two extrema of global and nearest-neighbor dispersal. In all cases, the tran- 
sition to coexistence occurs for values of the immigration rate well described 
by the neutral prediction, ud^N) 2/(2 -|- A^). Unfortunately, we could not 
systematically explore larger values of 7 or A^ as it requires huge statistics for 
distinguishing between the mixed and coexistence regimes. In fact at increas- 
ing 7 and/or A^ the distribution becomes strongly peaked on A^/2 with very 
low probabilities on the tails, so that the region where the differences between 
mixed and pure coexistence manifest is inaccessible (see Figj5]). 

Summarizing, the main effect of increasing 7 (from top to bottom in Figure H]) 
is to stabilize the coexistence by increasing the portion of the N — u plane 
where the mixed regime is realized. An equivalent analysis could in principle 
be performed by plotting the transition lines in the N — j and/or in the z/ — 7 
plane. However, the difficulty in sampling the tails of the distribution strongly 
limits the range of values of 7 and A^ which can be explored, so that these 
plots would not add much information to Figure S]). 

Neutral vs Non-neutral: coexistence & invasibility 

We now directly compare neutral and non-neutral dynamics by exploring 
whether the former (7 = 0) can reproduce/mimic the latter, e.g., in terms 
of generating similar patterns of coexistence such as the distribution P{Na)- 
As both 7 and u promote coexistence, it is natural to expect that the lack of 
habitat preference in the neutral model could be compensated by a larger im- 
migration rate u to reproduce similar distributions. However, as we will see, 
the different stabilizing effects of immigration and habitat preference have 
interesting dynamical consequences. 

In particular, we consider two attempts to reproduce non-neutral distributions 
with neutral ones for two different community sizes A^ = 100 (Fig. and 



11 



= 400 (Fig. |5hi). In both cases, we proceed as follows. We fixed the habitat 
preference intensity, 7 = 1 in this specific example. Then we simulated the 
neutral version of the model (7 = 0) and varied v until the distribution looked 
as similar as possible to that obtained with habitat preference. In other words, 
we searched for the value of the immigration rate v that compensated best 
the absence of habitat preference. As expected, this compensation is obtained 
by using a larger value of v. 

For = 100, the agreement between the curves is very good in the central 
peak, meaning that the differences between the two models can be appreciated 
only when one of the two species is dominating. For = 400 the distributions 
are almost indistinguishable. The fact that the non-neutral system of Fig |5^ 
is in the mixed phase is clear by looking at the distribution tails; however this 
feature could be difficult to detect in an experimental time series. In Fig. I^b, 
the tails have such a low probability that are essentially inaccessible even in 
long simulations, so that we cannot distinguish between mixed and coexistence 
phases. 

For large A^, both in the coexistence and mixed phase, the distributions around 
the central mode are well fitted by a Gaussian. In the neutral model, the 
limiting Gaussian distribution can be analytically derived (see Appendix |B|) . 
showing that the variance in this case is proportional to A^. Numerical sim- 
ulations (not shown) suggest that variance in the non-neutral version of the 
model scales in the same way with A^, with a prefactor rapidly decreasing at 
increasing 7 refiecting the presence of niches. 

Even if the difference between the models is essentially undetectable when 
looking at the distributions, it can play an important role when studying 
invasibility properties, i.e. by considering a situation in which one of the species 
is absent at the initial time (which means to directly test the tails of the 
distribution). As an example, for the same parameters of Fig. |5]we prepared 
the initial state with only species B present and we computed the average 
time it takes to reach perfect coexistence (A^^ = = N/2) for the first 
time. When A^ = 100 (a), this time is roughly equal to 46 generation for the 



Fig. 5. Non-neutral vs neutral distribution of population P{Na)- (a) for N = 100 
with 7 = 1, = 0.005 (solid) and 7 = 0, 1/ = 0.12 (symbols): (b) for = 400 with 
7 = 1, 1/ = 0.001 (solid) and j = 0, v = 0.123 (symbols). In both cases, the y axis 
is in logarithmic scale. 
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model with habitat preference and 17 generation for the neutral one, while for 
= 400 (b) the difference is even more pronounced (61 vs 21 generations, 
respectively). 



Habitat preference statistics 



We now investigate the role of preference sites in achieving coexistence and 
determining the spatial distribution of species. The natural quantity to look at 
is the average fraction of individuals of both species occupying their preference 
sites, Pocc = {NAa + NBb)/N . This quantity is a suitable measure of the non- 
neutrality of the system: it is equal to 1/2 in the neutral individuals 
have no reason to prefer one of the two site types. The average occupation 
increases with 7 eventually reaching 1 when the distribution of sites strongly 
determines the spatial distribution of the two species. 

The average occupation Pocc as a function of 7 is shown in Fig. [6] . For illus- 
trative purposes, here, we fixed = 100 and varied v. As one can see the 
curves are weakly dependent on the choice of z/, especially when 7 is large. 
The dependence on A^ is not shown but it is also rather weak. We remark that 
in the nearest neighbor case it is crucial to average over different realizations 
of the dynamics and of the distribution of advantageous sites: different spatial 
arrangements of the sites may be harder (or easier) to occupy for one of the 
species. 

When the two species are exactly equipopulated and with global dispersal, 
one can easily derive from ([2]) that the average occupation must be equal to 
(1 + 7)7(2 + 7). This prediction is shown for comparison in Fig. |6^. Significant 
deviation are present for small 7 and v values, for which pure coexistence is 
not achieved and species are not equipopulated, i.e. the system is in the mon- 
odominance or mixed phase. However, it must be noticed that with nearest- 
neighbor dispersal (Fig. [6]d) even for large 7 the prediction (1 + 7)/(2 + 7) 
is never realized due to dispersal limitation, which can prevent species from 
colonizing advantageous sites. In this case the average occupation does not 
reach 1 even for very large 7. The presence of patchiness in the distribution 
of preference sites is expected to enhance this effect: one could encounter sit- 
uations in which a whole patch is unreachable because it is surrounded by 
patches being advantageous to the other species. On the other hand, a longer 
dispersal range could compensate for the presence of patches. In other words, 
the possibility for the occupation to reach one for large 7 will depend on the 
ability of species to reach all the favorab le patches, which is known to be a 



crucial feature in fragmented landscapes (IKeitt et al.l . 119971 ). 



We conclude this section by discussing whether the results of Fig. [6] are de- 
termined by the particular choice of the ecological advantage. To check the 
robustness of our results, we simulated the variant of the model in which the 
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Nearest Neighbor Dispersal 




Fig. 6. Average preference occupation Pqcc as a function of 7 for N = 100 and 
several values of as in labels, (a) Results obtained with global dispersal, (b) near- 
est-neighbor dispersal. Data have been averages over 10^ generations. Data obtained 
with nearest-neighbor dispersal are also averaged over 200 independent realization 
of the initial distribution of preference sites. The continuous curve represents the 
simple prediction Pqcc = (1 + 7)/(2 + 7) (see text). 

advantage gives rise to a lower mortality in the preference sites, instead of a 
larger colonization probability, obtaining very similar curves to those of Fig [61 



4 Discussion 



We have introduced an individual-based model of competition between two 
species in the presence of habitat preference. As is typical in stochastic mod- 
els of community assembly, the system drifts towards the extinction of one of 
the two s pecies when immigration froni the metacommunity is not ta ken into 



account (IMacArthur and Wilsonl . 119671 : iLoreau and Mouquetl . Il999l ). Never- 



theless, even without immigration, habitat preference has a strong stabilizing 
effect on the coexistence in large communities, leading to a dramatic increase 
of the average time fo r the extinction of a species, consistently with studies of 
niche modelsjTilmanl, I2OO4J : lAdler et all . l2007l:lLin et al.l.l2009jK In particular 
Lin et al. fl2009h. extending previous results (jZhang and Linl . 119971 : lYu et al.l 



19981 : IWrightl . |2002| ) . introduced a model where both mortality and fecundity 



were varied among the species in such a way that the averaged life-history 
fitness of several species is the same, so to fulfill the weaker requirement of 
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the equivalence of average fitness of the neutral theory (IHubbelll . l2006l ) . The 
presence of birth/fecundity tradeoffs (globally equalizing the performances of 
different species) has a stabilizing effect, without affecting the ability of pre- 
dicting biodiversity patterns. 

The main message of the present study is that neutral and non-neutral coexis- 
tence may yield almost identical patterns such as the population distribution 
(see, e.g., FigjSjo) while dynamic pro perties may be sign ificantly different, as 



revealed by invasibility experiments (jPaleo et al.l . |2009[ ) . In this perspective. 



the advantage of the proposed model is to allow for quantitatively compar- 
ing niche and neutral competition in a simple setting and thus for providing 
an insight on the origin of such differences. It is thus worth discussing the 
predictions of the model for invasibility experiments. 

The simplest setting is that of an isolated ecosystem (i.e. without migratory 
flux), in which a new species with a small population is introduced. This 
corresponds to Figj2]D, where the probability of the invading species to take 
over grows with 7. The ecological interpretation is that niche-based assem- 
blies are easier to invade than neutral ones if the niches are not saturated, 
simply because the ava ilability of a free niche facilitates the invasion (see 
Melbourne et al. fl2007t ) for a recent review on the role of heterogeneity on 



invasibility properties) . 

Another possibility is to consider the invasion of an ecosystem initially pop- 
ulated by a dominant species in the presence of a regular flux of immigrants 
species from the metacommunity. In this case, if we compare neutral and non- 
neutral dynamics at equal immigration rate there is hardly any difference with 
respect to what is observed in the absence of immigration. 

Conversely, the response may be quite different if the comparison is made 
at equal realized distributions (see FigJS]). The neutral model realizes pat- 
terns of coexistence similar to non-neutral ones by compensating with a larger 
immigration pressure. Under these circumstances, if some catastrophic event 
leads to the extinction of one of the species, the non-neutral dynamics takes 
longer than the neutral one to recover the equipopulated state. Therefore, 
non-neutral coexistence is more robust, but it is also harder to achieve. The 
presence of two habitats makes both the two species coexist and one species 
hard to invade; in other words, niche-based coexistence is history-dependent. 
An important practical consequence is that inferring immigration rates from 
the fit of the distributions with neutral models could result in overestimates, 
even in situations in which the quality of the fit is very good. 

The above result agrees w ith the observation that large scale, specie s-rich 



patches are easier to invade ( jRobinson et al.l . ll995uStohlgren et all . 120031 ) . con 



trarily to the prediction of theories of competitive exclusion via niche partition- 
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ing and data from small scale observations (ITilmanl . Il997h. These contrastin g 



evidences constitute the so-called "invasion paradox" (IFridley et al.l . 120071 ) . 



Neutral theory predicts that the chances of establishment of an alien species 



depen d on the frequency of introduction of new individuals only (jPaleo et al. 
20091 ). On the other hand, in a non-neutral community this chance is strongly 



influenced by habitat diversity and availability of free niches. Our results sug- 
gest that, when comparing the invasibility properties of ecosystems, it is crucial 
to establish in which measure the diversity is maintained by the intensity of 
the immigration flux and/or by niche availability. The two-species framework 
considered here the stationary solutions have not t he same complexity as those 
of multi-speci e s mod els. However, previous studies (IChave et al.l . l2002uMcGilll . 
20031 : iTilmanl . l2004j ) also pointed out that species abundance distributions of 
neutral and niche-based multi-species models are very similar, provided suit- 
able parameters are tuned and suggesting that at least the qualitative features 
of invasibility dynamics here identified should hold also in the multi-species 
case. We recall that in multi-species models speciation, as well as immigra- 
tion, contributes to the introduction of new species, especially when large 
systems are considered. Therefore, the problem of inferring these rates from 
observed patterns becomes even more dramatic, since it is hard to compare 
these estimate with reliable measures of speciation rates. This could ex plain 
the high spec i ation r ates predicted by neutra l theor y (see iRicklefsl (120031 ) and 
also lHubbeli f l2003h : [Haeeeman and Etiennd tOO^ ). 



Finally, the occupation probability provides a measure of the departure from 
the static equilibrium case in which each species is confined to its preference 
sites. It is known that models correlating species distributions with environ- 
mental and climatic features cannot ful l y dete rmine the observed geographi- 
cal range of species ( IS Venning and Skovl . |2004| ). In fact, a complete predictive 
power of species- distribution models would r equire the existence of a static 
equilibrium state ( IGuisan and Thuillerl . 120051 ) . an assumption that is directly 
challenged by the neutral theory and partially violated by niche-neutral mod- 
els like the one presented in this paper. The departure from the equilibrium 
case is enhanced by dispersal limitation, as shown by the comparison between 
the nearest neighbor case and the global dispersal of Fig. [61 Indeed, this ob- 
servable is the one showing the most significant dependence on the choice of 
the dispersal among those considered in this work. This difference will become 
more dramatic if one adopts the more realistic choice of a correlated environ- 
ment, i.e. a distribution of preference sites being correla ted in space . As i t 



has been studied for the case of fragmented landscapes ( iKeitt et al.l . 119971 ) 



the comparison of dispersal ranges and characteristic size of patches deter- 
mines whether species will be able to reach all their preference sites or will be 
arrested by dispersal limitation. 



Summarizing, we studied the effect of habitat preference on the stochastic 
competition between two species. In the absence of immigration, demographic 
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stochasticity eventually leads to the extinction of one of the two species. How- 
ever, habitat preference leads to a dramatic increase of the extinction time. 
The drift to extinction is arrested by introducing an immigration rate. When 
the latter is large enough, the model generates a coexistence state in which 
both species occupy a significant fraction of the system. By compensating 
with a larger immigration rate the loss of specificity, the neutral case (7 = 0) 
can generate distributions which are essentially indistinguishable from those 
obtained in the presence of habitat preference (7 > 0). In spite of the pat- 
tern similarity, the dynamics of the two cases is very different. When only one 
species is present at the beginning of the simulation, the neutral variant of the 
model is much easier to invade than the non-neutral one. This latter result 
reproduces in a simplified framework what is observed in species-rich com- 
munity models: neutral and non-neutral models reproduce similar patterns of 
species abundances, here of population distribution, but can be distinguished 
by looking at dynamical properties. 
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A Condition for coexistence in the neutral model with global dis- 
persal 



For 7 = 0, i.e. in the absence of habitat preference t he model wit h global 
dispersal reduces to the Moran model with mutation ( iMoranl . 119581 ). In this 
case the transition rates can be expressed in terms of a single population Na- 
In particular, equations ([2]) reduce to the birth and death rates of population 
Na, i.e. W Af^^/v^+] = W ^(Na) which, incorporating the effect of immigration. 



read flKarlin et al.l . Il962[ ): 



W+{Na 



N-Na 
N 



;i-^)^ + f 



Na 
N 



(Al) 



N-Na 
N 



+ 



These rates can be used to calculate the equilibrium distribution P{Na) 
through the detailed-balance relation 



PiNA) 



W-iNA + l) 



(A.2) 
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which must hold at stationarity since the process is one dimensional (IGardiner 



20041 ) ■ The above expression can be used to determine the critical value of the 
immigration rate for the transition from monodominance to coexistence. It 
is enough to determine the value of v for which P{N a) passes from being a 
decreasing function of A^i (we are assuming A^^ < A^/2) to an increasing one. 
Coexistence is thus obtained when P{Na + 1) > P{Na), which using flA.ll) 
and after some algebra can be recast as 

[(2 + iV)z/ - 2] (A^ - 2A^A - 1) > . (A.3) 

Therefore, whenever i> > 2/{2 + N) the distribution is increasing up to Na = 
{N — l)/2 then decreasing, i.e. there is coexistence. When i/ < 2/(2 + A^) the 
opposite happens and there is monodominance. We conclude noticing that for 
u = 2/{2+N), i.e. at the the transition between the two classes of distributions, 
the distribution becomes uniform, i.e. P{Na) = l/N. This is true only for the 
neutral model with global dispersal. 



B Gaussian limit of the neutral model 



In this Appendix we show that the stationary solution approaches a Gaussian 
in the large A^ limit, at least in the neutral, global dispersal case. We perform 
the calculation in the diffusion approximation ; the procedure is s imilar to the 
classic results for the Moran model, see, e.g. (IKarlin et al.l . Il962l ). apart from 
small differences in how the immigration (mutation) mechanism is defined. 
From the rates flA.l|) the diffusion approximation yields the Fokker-Planck 
equation: 

dtp{x,t) = -drc{a{x)p{x,t)) + -dl{b'^{x)p{x,t)) 



where x = Na/N and 

u{l-2x) 
a{x) = 



2x(l - x)(l - z/) + z/ 
N ■ 



The stationary distribution is easily found to be 

Nu 1 

Ps(x) oc [2a;(l - a;)(l - z/) + z/]2(i— ) 2. 



(B.l) 



(B.2) 



(B.3) 



When A^ is large, an expansion around the maximum xq = 1/2 leads to the 
Gaussian distribution: 



Pst{x) 



1 



V2 



exp 



{x - Xof 



2a2 



(B.4) 



with 0"^ = (1 + i')/[A{N + l)z/ — 4]. This shows that the relative fluctuations 
around the maximum have the expected characteristic size 1/\/N. 
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